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With the puncture method for black hole simulations, the second infinity of a wormhole geometry 
is compactified to a single "puncture point" on the computational grid. The region surrounding the 
puncture quickly evolves to a trumpet geometry. The computational grid covers only a portion of the 
trumpet throat. It ends at a boundary whose location depends on resolution. This raises the possi- 
bility that perturbations in the trumpet geometry could propagate down the trumpet throat, reflect 
from the puncture boundary, and return to the black hole exterior with a resolution-dependent 
time delay. Such pathological behavior is not observed. This is explained by the observation that 
some perturbative modes propagate in the conformal geometry, others propagate in the physical 
geometry. The puncture boundary exists only in the physical geometry. The modes that propagate 
in the physical geometry are always directed away from the computational domain at the puncture 
boundary. The finite difference stencils ensure that these modes are advected through the bound- 
ary with no coupling to the modes that propagate in the conformal geometry. These results are 
supported by numerical experiments with a code that evolves spherically symmetric gravitational 
fields with standard Cartesian finite difference stencils. The code uses the Baumgarte-Shapiro- 
Shibata-Nakamura formulation of Einstein's equations with 1+log slicing and gamma-driver shift 
conditions. 



I. INTRODUCTION 

The puncture method [T, ^ for black hole simulations 
works remarkably well. Initially, each black hole is repre- 
sented as a wormhole. The second infinity of each worm- 
hole is compactified to a point on the computational grid 
called a puncture. The spatial geometry is evolved using 
finite differencing with the Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) formulation of Einstein's evolution 
equations [31 S] and the standard gauge conditions con- 
sisting of 1-t-log slicing and gamma-driver shift [6J. 




FIG. 1: Penrose diagram showing the relationship between 
the initial and evolved wormhole slices and a trumpet slice. 
The heavy dots represent the distribution of numerical grid 
points. This figure is a sketch based on the results from 

Refs. [2i[g. 

The evolution of a single, spherically symmetric punc- 
ture black hole is well understood from a geometrical 
point of view [7J [51 [HI HO] • Figure [l] shows the Penrose 
diagram for a Schwarzschild black hole. The initial worm- 
hole slice stretches between the two spacelike infinities. 
The puncture at r = coincides with the left spacelike 
infinity. Heavy dots mark the locations of points in the 



computational grid. The puncture itself is not a point 
in the computational grid, since the gravitational field 
diverges there. (Alternatively, one can apply a regular- 
ization scheme to keep the fields finite at the puncture 

[Hill [13].) 

The spatial metric for the initial wormhole slice is writ- 
ten as ds^ = ij^^^dr"^ + r^dil'^), where dil^ is the metric 
on the unit sphere. The conformal factor ip diverges like 
ijj ^ 1/r at the puncture. With l-l-log slicing and an 
initially constant lapse, the initial wormhole slice evolves 
to a slice that is left -right symmetric in the Penrose dia- 
gram. The gamma-driver shift condition, which is built 
from the conformal metric, breaks the left-right symme- 
try. This allows the gamma-driver shift to drive the grid 
points away from the left spatial infinity, from left to 
right in the Penrose diagram. 

The computational grid only covers a portion of the 
evolved wormhole slice. That portion is indistinguishable 
from a portion of a "trumpet slice" of Schwarzschild. A 
trumpet slice is a stationary l-l-log slice that asymptoti- 
cally approaches the left future infinity The confor- 
mal factor for a trumpet slice diverges like -0^1/ ^/r. 
Near the puncture r = the trumpet geometry consists 
of an infinite throat with topology R x and constant 
cross-sectional area. 

Only a finite portion of the trumpet geometry is cov- 
ered by the computational grid. The grid ends at the set 
of grid points that are closest to r = 0. The innermost 
layer of grid points can be viewed as a discrete boundary 
that divides the throat into interior and exterior regions. 
As discussed in Ref . f7| , the innermost layer of grid points 
constitute a natural excision boundary. The central issue 
addressed in this paper is how we should view this "punc- 
ture boundary". There is, of course, an outer boundary 
to the computational domain as well. The numerical rel- 
ativity community has spent, and continues to spend, a 
good deal of effort in designing boundary conditions for 
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the outer boundary. (See, for example, Ref. [H].) Should 
we also worry about boundary conditions at the puncture 
boundary? 

One might argue that the puncture boundary cannot 
affect the exterior evolution of the black hole since it lies 
entirely in the black hole interior. The situation is not 
that simple. We know from the analysis of hyperbolicity 
for BSSN with standard gauge conditions that there are 
a number of modes that travel with superluminal speeds 
[T51 [T^ . These modes represent the gauge freedom in 
the system, the freedom to change the slicing and spatial 
coordinates. In particular, these modes do not violate 
the constraints [13]. Nevertheless, it is possible that a 
gauge wave could travel into the black hole, reflect from 
the artificial puncture boundary, and propagate back to 
the black hole exterior. This raises a number of ques- 
tions. What are the boundary conditions at the puncture 
boundary? Are they determined by the finite difference 
stencil? How do the boundary conditions affect the re- 
flected waves? 

A more alarming issue is how the grid resolution might 
affect waves that reflect from the puncture boundary. As 
the computational grid is refined, grid points are added 
closer to the puncture point r = 0. The innermost layer 
of grid points is moved farther into the trumpet throat, 
so the distance is increased between the puncture bound- 
ary and any finite location, say, the black hole horizon. 
Of course, when we refine the grid the coordinate loca- 
tion of the puncture boundary changes by an infinitesimal 
amount proportional to the grid spacing. However, the 
proper distance from the horizon to the puncture bound- 
ary changes by a finite amount because the geometry 
diverges at r = 0. This raises the possibility that wave 
refiections from the puncture boundary will be delayed 
as the grid resolution is increased. 

In Sec. II I present simulations of a scalar field on a 
cylinder R x S"^ with puncture compactification. Spher- 
ically symmetric wave pulses are sent down the cylin- 
der and allowed to reflect from the puncture boundary. 
This system exhibits the pathological behavior described 
above. Speciflcally, the reflected waves show a clear 
resolution-dependent time delay. Moreover, we find that 
the form of the reflected wave is affected by the choice of 
flnite difference stencil at the puncture boundary. 

It appears that such pathological behavior does not 
occur in black hole simulations that use the puncture 
method. This paper is devoted to explaining why this 
is so. The picture that emerges can be summarized as 
follows. The initial wormhole data evolves very rapidly 
into a trumpet configuration. For the trumpet geome- 
try, only two of the characteristic fields are outgoing at 
the puncture boundary.^ These modes have sufficiently 



large positive speeds to allow perturbations to propagate 
from the puncture boundary to the black hole exterior. 
In effect, these modes propagate in the conformal geom- 
etry, not the physical geometry. As such, they do not 
sense the movement of the puncture boundary when the 
grid resolution is changed. The remaining characteris- 
tic fields have negative coordinate speeds at the punc- 
ture boundary. If a perturbation occurs in one of these 
modes, the time it takes to reach the boundary will de- 
pend on resolution. This raises the possibility that one 
of these incoming modes will be coupled to an outgoing 
superluminal mode at the puncture boundary and give 
rise to a resolution-dependent reflection. Numerical sim- 
ulations indicate that there are no such couplings. This 
can be understood by considering the finite differencing 
scheme as it appears in both the physical and conformal 
geometries. For the outgoing superluminal modes, the 
finite difference stencil is equivalent to a typical stencil 
that would be used to evolve smooth fields at the origin 
of a smooth geometry with topology R^. For the modes 
that propagate in the physical geometry, with topology 
R X S^, the stencil is one sided at the puncture bound- 
ary. These ingoing modes are simply advected off the 
grid with no coupling to the outgoing modes or to each 
other. 

Section III begins with a review of the covariant for- 
mulation of BSSN with standard gauge conditions [T7] . 
I then present a detailed derivation of the characteristic 
fields and speeds for this system. The analysis uses the 
"frozen coefficients" approximation, defined by a small 
amplitude, high frequency limit for perturbations on a 
background solution. In this way the characteristic fields 
(or "modes" ) are identified as perturbations that can be 
realized numerically. I also define a "gauge system" that 
has the same superluminal characteristics as the full sys- 
tem of BSSN with l-l-log slicing and gamma-driver shift. 
The gauge system is linear. Useful insights into the full 
system can be gained by examining the simple gauge sys- 
tem. 

In Sec. IV I show graphs of the characteristic curves for 
BSSN with the standard gauge, as an initial wormhole 
geometry evolves to a trumpet. Here we see explicitly 
which modes can carry information from the puncture 
boundary to the black hole exterior. Some of the modes 
are initially outgoing at the puncture boundary, but then 
quickly change to incoming. In principle, these modes 
should be fixed by boundary conditions while they are 
outgoing. In practice this does not seem to be a prob- 
lem, perhaps because these modes spend such a short 
amount of time with their characteristics outgoing at the 
puncture boundary. 

In Sec. V I describe the results of numerical experi- 



Throughout this paper I use the terms "outgoing" and "out- 
ward" to mean "in the positive radial direction; away from the 
puncture boundary". This differs from the common definition 



of "outgoing" as "from the interior to the exterior of the com- 
putational domain" . Likewise, terms such as "incoming" and 
"inward" will mean "in the negative radial direction; toward the 
puncture boundary" . 
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ments in which the characteristic fields are evolved as 
perturbations on a single, stationary trumpet geome- 
try. These experiments are used to probe the puncture 
boundary. The perturbations consist of a simple wave 
pulse of compact support in one of the incoming modes. 
The system is evolved and reflections in the outgoing 
modes are examined. This technique is used to identify 
any resolution-dependent time delay that might appear 
in the reflected waves. No such delay is found. 

Section VI contains a summary and discussion of the 
conclusions that can be drawn from the numerical exper- 
iments. 

The simulations presented in this paper use a cartoon- 
type code HSI based on the covariant version of BSSN de- 
scribed in Ref. [IT] . The cartoon code is described in Ap- 
pendix A. The code is designed to evolve spherically sym- 
metric gravitational fields using standard Cartesian coor- 
dinate finite difference stencils. In Cartesian coordinates, 
covariant BSSN is completely equivalent to the standard 
BSSN system. With the covariant formulation we can 
define transformations between Cartesian and spherical 
coordinates in a meaningful way. 

The cartoon code is designed to be third-order conver- 
gent. In Appendix A I present the results of a two-point 
convergence test with the Hamiltonian constraint and a 
three-point convergence test with the conformal factor. 
These tests confirm that the code is third-order conver- 
gent everywhere in the computational domain. It is often 
stated that puncture evolution codes do not converge at 
points near the puncture. Strictly speaking, this claim is 
incorrect. Although there are large finite differencing er- 
rors near the puncture, a properly constructed code will 
be convergent everywhere except at the puncture itself, 
where the conformal factor is not defined. 

Appendix B contains a detailed discussion of the two 
techniques that I use to generate the stationary trum- 
pet geometry for the numerical tests. I discuss how the 
numerical data is corrected to account for the fact that 
neither of these techniques yields a geometry that is pre- 
cisely stationary. 



II. PUNCTURE EVOLUTION OF A SCALAR 
FIELD 

Consider a massless scalar field (f> propagating on a 
three-dimensional cylinder with topology R x S^. The 
metric is ds"^ = d^'^ + dd^+sm^ 9 dcjp where — oo < C < oo 
and 9, (p are spherical coordinates on 5^. Now compactify 
the second infinity ( = — oo with the coordinate trans- 
formation 

This transformation maps the three-dimensional cylinder 
to with a puncture at the origin r = 0. The spatial 



metric gab becomes 

ds^ = (^1 + ^) dr^ + d9^ + sin^ 9 dcj)^ . (2) 

This metric agrees with the trumpet metric in the limit 
r 0. In particular, both describe a three-dimensional 
cylinder with constant cross-sectional area and proper 
length that diverges like l/r. 

The scalar field equation with unit lapse and vanishing 
shift is 

= n , (3a) 
dtU = g'^'-i^aA* . (3b) 

Here, Da is the covariant derivative built from the metric 

©■ . . . 

Numerical simulations are carried out using the Carte- 
sian coordinates defined by r = x"^ + y'^ + z^, cos 9 = 
z/r, and tan^ = y/x. The puncture resides at the co- 
ordinate origin. The numerical grid points lie at the in- 
tersections of the coordinate lines x/h = 1/2,3/2,...; 
y/h = 1/2,3/2,...; and z/h = 1/2,3/2,...; where 
h is the grid spacing. We can visualize this grid by 
drawing the coordinate lines x/h = 1/2,3/2,... and 
y/h = 1/2,3/2,... in the equatorial plane defined by 
z = 0. Figure [2] shows these coordinate lines plotted 
on a graph of C versus 0. The left and right edges 
are identified, so the figure represents an infinitely long 
two-dimensional cylinder with circumference 27r. In this 




FIG. 2: Coordinates on the compactified cylinder. The left 
and right edges of this figure are periodically identified. The 
grid points lie at the intersections of the coordinate lines. 

two-dimensional figure, the discrete puncture boundary 
consists of the four points {x,y) — {±h/2,±h/2) that 
coincide with the four intersections near the bottom of 
the figure. In the actual three-dimensional simulations, 
the puncture boundary consists of the eight grid points 
{x,y,z) = (±/i/2, ±/i/2, ±ft,/2) that are closest to the 
origin. 
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Figures [3] through [6] show a sequence of images taken 
from simulations at three different resolutions, low (h = 
1/25), medium {h = 1/50), and high {h = 1/100). The 
initial data consists of a spherical Gaussian pulse 



^ g-2(r-5)^ ^ 



n(0) = -4(r-5)e 



(4a) 
(4b) 



that is traveling toward the puncture r = 0. At time 
t = A, the pulse is still traveling toward the origin and 
the wave forms obtained from the three simulations coin- 
cide. The wave pulse hits the puncture boundary around 
time t — 8. The figures at times t = 12 and t = 16 
show the reflected pulse propagating away from the ori- 
gin. There is a clear dependence on resolution, with the 
low resolution pulse being reflected most quickly. The 
reflection of the medium resolution pulse is delayed rela- 
tive to the low resolution case, and the reflection of the 
high resolution pulse is delayed even further. 
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FIG. 4: Scalar field propagating on a punctured cylinder at 
time t — 8. The wave is interacting with the puncture bound- 
ary. 
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FIG. 3: Scalar field propagating on a punctured cylinder at 
time t — 4. Data for three resolutions are shown. The three 
curves overlap at this time. 

The resolution dependence of the reflected pulse is pre- 
cisely what we expect, based on the following analysis. 
Scalar flelds propagate with proper speed di/dt — ±1, 
where d£ is proper length. Then the coordinate speeds 
in the radial direction are dr/dt — zLl/y/g^. With the 
metric (jij), the wave speeds are Vin = — r/vT+T^ for in- 
coming waves and Vout — r/^/T+V^ for outgoing waves. 
The time required for a pulse to propagate from some 
large radius rg to the innermost grid point at r ^ h/2 
and return to rg is 



h/2 



1 



dr 



1 

h/2 '"out 



dr 



const — 2 \og{h) . 



(5) 



The constant in this expression is independent of grid 
spacing h. It follows that the time difference between 



FIG. 

time 
ture. 



5: Scalar field propagating on a punctured cylinder at 
t = 12. The reflected waves are emerging from the punc- 



resolutions h and h/2 is Th/2 - T^. = 21og2 « 1.37. We 
can compare this with our numerical results. At time t — 
20, the reflected Gaussian pulses are peaked at r « 8.35, 
6.97, and 5.61 for low, medium, and high resolutions, 
respectively. The difference between low and medium 
is Ar Ri 1.38, whereas the difference between medium 
and high is Ar « 1.36. At these radial distances the 
wave pulses travel with coordinate speeds very close to 
unity. Thus, the spatial separations of 1.38 and 1.36 are 
remarkably close to what we would expect given a time 
lag of 1.37. 

It is important to keep in mind that the reflection is 
artiflcial. The compactiflcation has introduced an inner, 
"puncture boundary." For the ideal problem of a scalar 
field propagating on an infinite cylinder, the wave pulse 
would continue to travel down the cylinder with no re- 
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FIG. 6: Scalar field propagating on a punctured cylinder at 
time t = 16. The refiected pulse shows a distinct time lag 
that increases with resolution. 

flection. We can approach this ideal behavior by increas- 
ing the resolution. When we increase the resolution the 
puncture boundary is pushed farther down the cylinder 
and the artificial reflection is postponed. In principle, the 
artiflcial reflection can be postponed indeflnitely. Unfor- 
tunately, this is not feasible in practice. For example, if 
we want the reflected pulse to be delayed by a time of, 
say, 20, then Eq. ([5| shows that the required resolution 
is /i « e-i° « 1/22000. 

Now consider how the flnite difference stencil affects 
the reflected pulse. The simulations presented above use 
standard centered fourth-order stencils. For example, for 
the first and second derivatives of $ with respect to the 
coordinate z, we have 

+8$fc+i - $fe+2) , (6a) 
{dl<^)k = ^(-$fc-2 + 16$fe^i-30$fc 

+16$fc+i - $fc+2) , (6b) 

where k labels grid points. These finite difference deriva- 
tives are fourth-order accurate. Alternatively, we can 
compute the derivatives of $ by 

(9.$)fc = ^(-3$fc_i-10$fc + 18$fc+i 

-6$fe+2 + $fc+3) , (7a) 

(9'$)fc = ^(ll$fe-i-20$fc + 6$fe+i 

+4$fe+2 - $fc+3) • (7b) 

This first derivative is fourth-order accurate; the second 
derivative is third-order accurate. 

Figure [7] shows the results of a simulation in which 
the finite difference stencils are altered for the eight grid 



points closest to the puncture. These are the eight grid 
points that form the discrete puncture boundary. Sten- 
cils of the form ([?]) were used for first and second deriva- 
tives along the x, y, and z coordinate lines. Cross deriva- 
tives such as dxdy^ were not changed. 
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FIG. 7: Reflections from the puncture boundary of a scalar 
field pulse at three resolutions. The solid curves use the stan- 
dard centered finite difference stencils. The dashed curves use 
alternative stencils for the innermost layer of grid points. 

The solid curves in Fig. [7] are a repeat of the curves 
from Fig. [6| The dashed curves are obtained with the al- 
tered stencils applied to the innermost grid points. There 
is a distinct difference between these curves. With the 
altered stencils, the reflected pulses are shifted in time 
and altered in shape. This suggests that for the punc- 
ture evolution of a scalar field, the boundary conditions 
at the puncture boundary are determined, or at least af- 
fected, by the finite difference stencil. 



III. PERTURBATION ANALYSIS FOR BSSN 
AND THE STANDARD GAUGE 

A. Covariant BSSN 

The BSSN equations with l-|-log slicing and gamma- 
driver shift conditions are written in covariant form in 
Ref. [17], based on earher work in Refs. [19l |20]. The 
BSSN variables ip, gab, K, and Aat are defined in terms 
of the physical metric g^b^^ and extrinsic curvature K^^^^ 

by 

= ^'"gab , (8a) 
^Phys ^ e^^Aa, + ga,K/3) . (8b) 

Let AF;^^ = FJJ^ - t'^^ where F^J^ are the Christoffcl sym- 
bols built from the conformal metric and F^^ arc the 
Christoffel symbols built from a background metric. Also 
use the notation AF" = .g'^^AF"^. The BSSN variables 
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include the "conformal connection vector" defined by The BSSN and standard gauge equations in covariant 

form are p/7j 

A" = Ar" . (9) 
In this paper I will assume that the background is flat. 



dt9ab 
dtAab 

dtK 



iS'^d.Aab + 2A,^adb)l3' ~ lAabD,/3^ ~ 2aA,,Ag + aA^bK 



+e" 



-4ip 



[-2aDaDbip + AaDaipDbip + 4:D(^aaDf,jip - DaDba + aUab] 



TF 



6 



P'^dcK + -K^ + aAabA"^ - e"^'^ {D^a + 2D''aDaf) , 
o 



-2A''-{6^d,a - 6a6^d,^ - aAF^J - -ag'^'dbK , 



dta = P'^Dca ~ 2aK , 



where 



Tiab = -l^g" DcDdQab + gc{aDb)A'' + g ''Ar5eAr(ab)c 



cd 
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Ar 



b)ed 



ArLAr 



ebd 



(lOa) 

(10b) 
(10c) 
(lOd) 

(lOe) 
(lOf) 
(lOg) 
(lOh) 

(11) 



In the equations above it is assumed that A = g'^^Aab = is enforced. The superscript TF denotes the trace-free 
part of the expression in square brackets, defined with respect to either the conformal metric or the physical metric. 
The term (9tA°)rhs in Eq- (lOi) must be replaced by the right-hand side of Eq. (10^). The operator Da denotes the 
covariant derivative with respect to the conformal metric gab- The operator ba denotes the covariant derivative with 
respect to the initial conformal metric, gab(O). The operator ba denotes the covariant derivative with respect to the 
(flat) background metric. 

If the initial conformal metric is flat and the coordinates are Cartesian, then the equations of motion (lOi-e) 
above are identical to the original BSSN equations [3, 4 . Likewise, Eq. (10') is the usual 1+log slicing condition 
and Eqs. ( 10 5,h) are the gamma-driver shift conditions. Note that the gamma-driver shift equations include all 
advection terms. With this choice, Eqs. (10 5,h) are equivalent to the single equation spatial gauge condition used by 
the Goddard group |21j . In the final section I discuss how the results are changed if one or more of the advection 
terms are dropped. 

For the BSSN formulation of Einstein's equations the constraints are defined by 



n = -if^-^,fcA'^'' + e-^'^(i?-8Z?VI?,¥>-8i^V) , 

Ma = 9'''bbAac~AabAT''-A\ATlb + (!>Ald,^~'^daK 
C" = A" - AF'^ , 



(12a) 

(12b) 
(12c) 



where R is the conformal Ricci scalar. Here, 7i = and A^q = are the usual Hamiltonian and momentum constraints. 
The constraint C = arises from the definition of the confor mal connection vector. 



B. Characteristic fields tives can be a analyzed using the pseudodifferential oper- 

ator method of Refs. [22l|23]. This technique was first ap- 

Hyperbolicity for systems of partial differential equa- 
tions with first order time and second order space deriva- 
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plied to BSSN with standard gauge conditions by Beyer 
and Sarbach |l5j . 

It will be useful to begin by adopting a condensed no- 
tation. The BSSN plus standard gauge equations (10) 
have the quasilinear form 

dtq = A''{q)daq + B{q)v + Ciq) , (13a) 
dtV = D"-\q)dad,q + E\q)daV + F{q,dq,v) {lih) 

where q represents the "coordinate variables" gab, ot, 
and /3° and v represents the "velocity variables" Aab, 
K, A", and B". The coefficients A", B, C, D"-'', and 
depend on the q's. The function F is a quadratic 
polynomial in the variables daq, v with coefficients that 
depend on the g's. 

Let q, V denote a solution of Eqs. (13) and consider 
perturbations q, v about this solution: 



q = q + q , 

V — V + V . 



(14a) 
(14b) 



The "frozen coefficients" approximation is defined by a 
small amplitude, high frequency limit as follows. The 
coefficients in Eqs. ( 13 1 and their spatial derivatives are 
assumed to be of order unity or smaller. The perturba- 
tion fields have small amplitudes proportional to e and 
large wave numbers proportional to k. To be precise, we 
let dq ^ V ^ 0{e) where d denotes a space or time deriva- 
tive. For higher order derivatives, ddq ~ 9w ~ 0{e ■ k). 
The coordinate perturbation q satisfies q ~ 0{e/k). In 
the limit as e — > and fc — > oo, the leading order (nonva- 



nishing) terms in Eqs. ( 13 1 are 



dtq = A-{q)daq + B{q)v , 

dtV = D'^\q)dadbq + E-{q)dai 



(15a) 
(15b) 



Equation (15i) is derived from the 0{e) terms in 
Eq. (13i), and Eq. (15d) is derived from the 0{e ■ k) 



terms in Eq. (13 d). Equations (15) define the principal 



part of the system (13). Note that the coefficients A", 



B, D"'^ , and are functions of the unperturbed coordi- 
nates q. Below, I will drop the hats from the unperturbed 
fields for notational simplicity. 

The hyperbolicity of the system is analyzed by set- 
ting the perturbation to a single Fourier mode with wave 
number ka'- 



x) = ve 



(16a) 
(16b) 



Here, |fc| = ^/fcoS"^^ and gab is the unperturbed confor- 
mal metric. The unit normal to the wave front is denoted 
Ua = ka/\k\. The indices on ka and Ua are raised and 
lowered by the unperturbed conformal metric. Note that 
Eq. imphes n°-daq = ge*'^*+*'=''^° . 

Now let /i = w/l/cj. The linear Eqs. (15) for the per- 



turbations become a set of linear algebraic equations for 
the Fourier coefficients: 



Hq — A'^Ua q + Bit , 

— D°-^nanbq + E'^naV 



(17a) 
(17b) 



The principal symbol V is the matrix defined by the 
right-hand side of these equations. 



V = 



B 



D'^^naUf, E^Ua 



(18) 



The system of differential equations ( 13 1 is strongly hy- 



perbolic if V possesses a complete set of eigenvectors with 
real eigenvalues. 

The characteristic fields for the system are obtained 
from the left eigenvectors of the principal symbol. Let 
(^, C) denote such an eigenvector; that is, (f , C)V — 
/i(^,C)- The characteristic field associated with this 
eigenvector is x = ^n'^daq + Cv. It satisfies dtX = fJ-n'^daX 
in the frozen coefficients approximation. 



We now spell out the results explicitly for the BSSN plus standard gauge equations. The principal parts of these 



equations are 



dtgab 




4- 2g,^adb)P' 


2 ~ ~ 
- i^gabdcP" - 2aAab , 




(19a) 


dtAab 


= f3'd,Aab 




-2adadb(p - ^a^ba - '^g^'^dcdagab + a9c(adb)A'' 


TF 
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(19b) 




= I3'd,>f + 


\dj'^ - 
6 


- \ak , 
6 




(19c) 


dtk 


= P'dck- 






(19d) 








- Ig^'dhdJ' - -ag-^'dkk , 




(19e) 


dto. 


= P'dca- 










(19f) 


dtP" 












(19g) 




= pedes'" - 


f g'^dhdj- - 


H Ig^'dhdJ' - '^ag'^'d.k . 




(19h) 



Recall that the hats have been dropped from the unperturbed solution. As above, we define = \/kag°'^ki, and 
Ua = fca/|fc|, so that Ua is normalized with respect to the unperturbed conformal metric: Uag'^^rib — 1. The principal 
symbol is defined by 



^J■gab 


= P"9ab 


-f 2n^aPk 


2 

) - o5a6/3" - 2aAab , 




(20a) 








^2ana 


nbf ~ naTiba ~ -gab + ""-(aAfc) 


TF 


(20b) 




= + 


-P" ~ ^ak , 
6^ 6 






(20c) 


fj.k 


= p"k- 










(20d) 




= P^A""- 


^P^+l 


n'^P" - 


4 
3 




(20e) 




= /3"a- 


2ak , 








(20f) 


tip- 


= P^P"- 










(20g) 




^ P^B"- 


o 


4 
o 




(20h) 



Here, the notation T" = T°-na is used for any tensor T'^ contracted with the normal covector Ua- 

Define e\ by nae\ = and e\gab^B ~ ^ab, where the upper case indices A, B . . . range over the values 1 and 2. 
Thus forms an orthonormal diad in the subspace orthogonal to Ua- Indices on Ua and are raised and lowered 
with the conformal metric and its inverse. For any tensor we define T„ = Tan"" = T°'na = and T4 = T^e^. 
The diad index A is raised and lowered with the identity tensor 6ab and its inverse. 



With this notation we can split the principal Eqs. (20 1 into scalar, vector, and trace-free tensor blocks. The scalar 
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block is 



1^9 Tin 


= P^9nn + Pn ' 




(21a) 






(21b) 






' 4 ^ 2^ 2 ^ ' 

-gOf'P - 3" - ^9nn + 2 "An + ^9AA 


(21c) 


flip 


6 6 


aK , 


(21d) 


Ilk 


= p"-k - e-'^'^a , 




(21e) 


/"An 


= /3"A„ + - 


4 


(21f) 




= P"a - 2ak , 




(21g) 




= /3"/3„ + '^-Bn , 




(21h) 




= /3"B„ + pn - 


4 - 


(21i) 



where 9/1/1 = cjAB^^^ ■ In deriving these equations I have used the fact that the perturbation Aab is trace- free. This 
implies Ann + Aaa = so that Aaa can be eliminated in favor of Ann- The vector block of the principal symbol is 



fJ'gnA 
IJ'AnA 

IiKa 

(J-^A 

hBa 



P'^gnA +I3a- 2aAnA 



= P^AnA + e-^^ 

= /3"Aa + h , 
= + \ba , 

= P^Ba + ^a . 



:9nA 



a 



A. 



(22a) 
(22b) 

(22c) 
(22d) 
(22e) 



Finally, the trace-free tensor block is 



^^9^B = P'^gAB - 2a^AB , (23a) 
A = (3-A% - ^e-^^g'/^ , (23b) 



where the trace-free part of a tensor Tab is defined by T^^ = Tab — {TcdS'^^)Sab/'^- 
The characteristic fields obtained from the scalar block are 



XI = - 8d„^ , (24a) 

X2 = A" - 8dn^ , (24b) 

X3 = dnfinn + 3„5AA , (24c) 

Xt = -4re-^^dna ± k , (24d) 

xf = T^Ann ±k + e-^v'A" + ^e-^-^dngAA - le-^''dn~gnn - 2e-^'^dnV , (24e) 

= ^ (1 - 2ae-^^)B" ±ak + ae-^^a„a T (1 - 2ae-^'^)9„^" . (24f) 
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These relations can be inverted for the field perturbations as long as (1 — 2ae ^'^) 7^ 0. In that case we have 

,+ 



dnQAA 
A 

k 

A" 

dnCt 



XI + 3X2 +3X3 ^{^a) (i_2ae-4^) 3 ^ ' ^ [l - 2ae-^^) 

4 2 



Xi 



7X2 



:X3 



1 

3-^" ' 3 

1. + -^ 1. - 
3(xr-X4)+ glXs 



(2a)^/^e--^M^t^ + \e^nxt + X,^) " iMP^ 



3(1- 2ae-4¥') 



■xg 

_ 1(2^)3/2^-2. (X4+ + X4) 1 (Xg- + X6) 

8 24^ ^ (l-2ae-4v) 12(1-206-^"^) 



(X4+ 



-X4 



-Xi+X2-:^(2a)^/'e-2^ 



=) 12 (1 - 2ae-^v) 

ixt + Xi) , 2 {xt + X.e) 



(1 - 2ae-4v') 3(1- 2ae~^f) 



2 -e'^(X4+ + Xr) , 

« (x|-xr) _ 1 ixt - Xe ) 
2(1- 2ae-4v') 2 (1 - 2ae-4v) ' 

x3/2„-2cp (X4+Xr) 2_(X61±X61_ 

3 ^ ^ (1 - 2ae-4¥') 3(1- 2ae-4'P) 



(25a) 
(25b) 
(25c) 
(25d) 
(25e) 
(25f) 

(25g) 
(25h) 

(25i) 



The characteristic fields from the vector and trace-free 
tensor blocks are 



X7 

,± 



Ba~^a 
2 



Xf = BaT ^^riPA 



Xg = Aa - dn9nA T 2e^'^ii„A 

Xio 

The inverse relations are 



(26a) 
(26b) 

(26c) 
(26d) 



dngnA = 


-X7 + '^{xt + Xs - Xg - Xg ) , 


(27a) 


A-nA — 


^e"^'^(Xg" - Xg ) , 


(27b) 


Aa = 


1. + 

-X7+ 2(X? + X8) > 


(27c) 


dnPA = 


V3 _ + 
^(Xs -Xl) , 


(27d) 


Ba = 


^(X^ + X8') , 


(27e) 




e^'^ixto - xFo) ' 


(27f) 


Atf - 
^AB — 


^(xto + Xw) ■ 


(27g) 



Note that modes X7, Xg , and Xg are vectors in the sub- 
space orthogonal to Ua', the index A is suppressed above. 
Likewise, the tensor indices AB and trace-free symbol tf 
have been suppressed on Xw- 

The perturbations in the lapse and shift (and auxil- 
iary field) are built from modes X4 and Xe ■ These can 



be viewed as gauge modes. They do not violate the con- 
straints: In the frozen coefficients approximation the per- 
turbations of the constraints are given by 



n = 



le-^''dnX2 - le-*^dnX3 



Mn 

Ma 
Ca 



3e-'^a„(x5+ + X^) , 
dn{x^ - X^) , 

e~^^a„(x!r -x^) , 



rX2 



^X3 + i.e^'^ixt 



X5 ) 



^(X^ + X9 ) 



(28a) 
(28b) 
(28c) 
(28d) 
(28e) 



These perturbations depend only on modes X2, X3j Xs^j 
Xg and their spatial derivatives. 



C. Characteristic speeds 



The phase angle for the Fourier mode (16) is tot + 
kax'^ = \k\{ij.t + Uax""). The surfaces of constant phase 
satisfy Uadx"" /dt = — /i. Let us choose the wave fronts to 
coincide with the surfaces r = const, where r is one of the 
spatial coordinates. Then the covector ka is proportional 
to 5^ and Ua = S^^/^/g^^ . The coordinate speed of the 
wave is dr/dt, or 



coordinate speed 



-MV.9 



(29) 
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TABLE I: Coordinate and proper speeds for the scalar modes of BSSN with 1+log slicing and gamma-driver shift. 





Xi 


X2 


X3 


xt 


xt 


xt 


coordinate speed 










-13'' ± ae-^-^^/Tr 


-13'' ±^/g^ 


near trumpet puncture 


-0.53 r 


-0.53 r 


-0.53r 


-0.53 r± 0.72 r^/^ 


-0.53r ±0.35r2 


-0.53r± 1 


near wormhole puncture 











±5.7 


±4.0 


±1 


proper speed 











±-s/2[a 


±1 


±e2'^/a 


near trumpet puncture 











±2.1/^ 


±1 


±2.9/r2 


near wormhole puncture 











±1.4 


±1 


±0.25/r^ 



This is the wave speed as seen by the "Lagrangian ob- 
servers" who move along the d/dt coordinate lines. 

For "Eulcrian observers" who arc at rest in the space- 
like hypersurfaces, the wave speed differs by the addi- 
tion of the shift. To be precise, recall that —(3°'dt is 
the change during time dt in the spatial coordinate lo- 
cation of an Eulerian observer. Then the coordinate 
distance that the wave front travels in coordinate time 
dt, as seen by an Eulerian observer, is the difference: 
dr - dt) = dr + (3''dt. 

The proper distance or time between two surfaces 
a and a + da is ds = da j \J ±dacrj°-''dba, where 7°'' 
are the contravariant components of a spatial or space- 
time metric. This result is derived by first construct- 
ing the unit normal to the a = const surfaces: n" = 
7"*96CT/-y±^J<7Y^^9dfT. The unit normal can be writ- 
ten as n" = dx°/ds, where s is proper distance or 
proper time. Now compute da/ds = daa{dx°'/ds) = 
\/ ±daa'y°'^dba and the desired result follows. If a is one 
of the coordinates, the proper distance or time becomes 
ds = da/y/±Y^. 

The above argument shows that for an Eulerian ob- 
server, the wave front travels a proper distance {dr + 

/5''a'0/v/.9phys in a proper time of d^/^^^^- Here, 
S'sptm is the contravaxiant tt component of the spacetime 
metric, related to the lapse function by ffsptm = — 1/a^. 
Putting this together we find that the proper speed of a 



given mode is {dr/dt + P'')/{a^/9^s), or 

proper speed = -^^-j={f3'' - liy/g^) . (30) 

Note, r can be any one of the spatial coordinates. We will 
often choose r to be the radial coordinate in a spherical 
coordinate system. 

The speeds of the various modes are listed in Tables 
I and II. The first line in each table lists the coordinate 
speed of the mode in general. The second line of each 
table shows the coordinate speed of a radial wave near 
a trumpet puncture. These speeds are obtained from 
numerical data, which yield a « 0.46 r, (3^ « 0.53 r, 
g-2(^ ^ Q 75 J. grr ^ j^gj^j. 7. = for a single, spher- 
ically symmetric trumpet with M = I. (The numerical 
coefficients should be accurate to within 10%.) The third 
line of each table lists the coordinate speed of a radial 
wave near a wormhole puncture. These speeds are found 
from the initial data for a single, spherically symmetric 
black hole with M = 1 in isotropic coordinates. They 
also assume the initial conditions a = 1 and /?" = for 
the lapse and shift. Then near a wormhole puncture, we 
have a w L ,,3'' = 0, .9^'' = 1, and e^^v ~ 4^2. The 
fourth line of each table lists the proper speed of each 
mode. The fifth and sixth lines show the proper speeds 
near trumpet and wormhole punctures. 



The coordinate grid moves with respect to the Eulerian observers. The coordinate speed of this motion is /?'', and 
the proper speed is /{ce^/gl^) = e^'^ / {ctV 9^^) ■ Near the puncture boundary of a trumpet slice, the proper 
speed is 1.5/r. Near a wormhole puncture boundary the proper speed is zero, assuming the shift vector vanishes. 

D. Gauge System tions: 





dtK 




(31a) 




dta 


= p^d^a - 2ak , 


(31b) 




dtp'' 


= p^d.'p'^ + ^ 


(31c) 


Many of the issues that require close examination re- 
side in the "gauge sector" . We define the gauge sector 
by the following system of linear partial differential equa- 




= P^d,B'' + g'"'bbDeP'' 


(31d) 
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TABLE II: Coordinate and proper speeds for the vector modes (x7, ^ind x^) and trace-free tensor modes (Xk)) '^^ BSSN 
with 1+log slicing and gamma-driver shift. 





X7 




xt 


Xu) 


coordinate speed 








-13'' ±ae-^'-<',Jg^ 


near trumpet puncture 


-0.53r 


-0.53r±0.87 


-0.53 r± 0.35 


-0.53 r± 0.35 


near wormhole puncture 





±0.87 


4.0 


4.0 


proper speed 





±^e2^/(2a) 


±1 


±1 


near trumpet puncture 





2.5/r2 


±1 


±1 


near wormhole puncture 





0.22/r2 


±1 


±1 



Here, the fields g""^ , (p, a and are frozen. The principal 
part of this system coincides with Eqs. (191,f,g,h). The 
characteristic fields are Xa ^^id ■ The inverse relations 
are written in Eqs. (25^,g,h,i). 



This system is simple enough to allow us to compute 
explicitly the time evolution of the characteristic fields, 
with the following results: 



= Lt(.Xi) , 



')drXt 



= Lq (X4,X6) 



(32a) 
(32b) 



that come from the nonlincarities in the equations. What 



Eqs. (32 1 reveal is that, already at the level of the linear 



system (31), there is coupHng between modes Xa ^^'^ 
Xg . In particular, note that modes Xa ^i^l tend to excite 
modes ■ This coupling plays an important role in the 
analysis of Section V. 



geometry evolves from wormhole to trumpet. In each 
figure the dashed curve is the black hole horizon. Figure 
pi shows the characteristic curves for xt ^^^^ begin at 
time i = at radii r = 0.025, 0.05, . . . 1.0. The modes xt ^ 
Xg , and Xio travel with the speed of light. For these 
modes Fig. |9] shows the characteristic curves that begin 
at time t = and radii r = 0.025, 0.05, . . . 1.0. Figures [TOl 
and [TT] show the characteristic curves for the outgoing, 
superluminal modes xt ^'^^ xt ^ respectively. In both of 
these figures, one set of curves begin at time t = and 
radii r = 0.025, 0.225, . . . 0.825. Another set of curves 
begin at radius r = 0.025 and times t = 0.2, 0.4, . . . 4.0. 



These are the usual advection-type equations plus 
sources. The sources arise from the fact that the char- 
acteristic speeds are not constants. The source functions 
are linear in the characteristic fields xj with coeffi- 
cients that depend on the fixed fields 5"'', a, (3^ and 
their spatial derivatives. The source functions are 
linear in the characteristic fields X4 Siiid xt with coeffi- 
cients that depend on the fixed fields 5°'', ip, a, and ^ 
their spatial derivatives. 

For the full BSSN plus standard gauge system, the 
time derivatives of x^ ^ind xt include the same terms 
displayed in Eqs. (32 1. They will have extra source terms 




FIG. 8: Characteristic curves for mode xt ■ The dashed curve 
is the black hole horizon. 



IV. CHARACTERISTIC CURVES 

Consider a single, spherically symmetric 
(Schwarzschild) black hole evolved with the punc- 
ture method. The initial data is a wormhole in 
Cartesian coordinates, so the BSSN variables are 
(p = ln(l + l/(2r)), gab ^ Sgb, K = 0, Aab = and 
A" = where r — \/ + y"^ + z^. The lapse, shift and 
auxiliary variable are given initial values a = 1, = 0, 
and B"^ — 0. The black hole mass M is set to unity. 

Figures |8] through [TT] show the characteristic curves 
for the various modes from time i = to t = 10, as the 



It is clear from these figures that perturbations in 
modes xt '"^^^ X^ propagate from the puncture 
boundary to the black hole exterior. Figure |8] shows that 
mode xX superluminal; however, perturbations 

within a radius of about r ~ 0.2 do not escape to the 
outside. This result was observed in Ref. |13j . 

Tables I and II show that modes X4 , xt ^ Xg"; and Xio 
have positive coordinate speeds near the wormhole punc- 
ture. This is difficult to see in Figures |8] and |9] What is 
clear from these figures is that the positive speeds do not 
last for long. Very quickly, as the geometry shifts from 
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FIG. 9: Characteristic curves for modes xt j xt ^^'^ Xi 



FIG. 11: Characteristic curves for mode xt ■ 





FIG. 10: Characteristic curves for mode xt 



FIG. 12: Time for characteristic speeds to become negative 
as a function of the radius. 



wormhole to trumpet, the speeds become negative near 



the puncture. Figure 12 is a graph of the time required 
for modes xt > xt ' XgT and Xio acquire a negative 
speed. The time is plotted as a function of r, which 
should be viewed as the puncture boundary radius. For 
example, with a grid spacing of h = 1/25 the puncture 
boundary has a radius of r « 1/50 = 0.02. Then Fig. 12 
shows that mode xt ^ positive speed at the punc- 
ture boundary up to time t sa 0.3. Beyond t w 0.3, the 
characteristic speed for xt negative. 

In principle, any characteristic mode that can prop- 
agate into the computational domain (in our terminol- 
ogy, any mode that is outgoing at the puncture) should 
be fixed by boundary conditions. Typically, a numerical 
code that fails to fix such a mode will be unstable. Thus, 
one would expect that boundary conditions are needed 
for modes xt ' xt ' xt j ^^id xto during the early stages of 
their evolution. Yet, the puncture evolution scheme does 
not obviously provide any such boundary conditions. In 



practice this does not appear to matter. It seems likely 
that these modes have positive speeds for such a short 
amount of time that instabilities do not have a chance to 
grow. The poor resolution near the puncture boundary 
might keep the instabilities' growth rates very low. 

Consider again a simulation with h = 1/25. In this 
case the modes xt ^ xt ^ Xg ' ^^'^ xto have positive speeds 
and are "unfixed" for a time of i « 0.3. If the Courant 
factor is, say 1/4, then the numerical code will execute 
about 30 time steps before the speeds change sign. This 
is perhaps too few time steps for an instability to grow 
to a significant level. 



V. PERTURBATIONS OF THE TRUMPET 
GEOMETRY 

In this section we carry out simulations of waves re- 
flecting off the puncture boundary of a trumpet black 
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hole geometry with M = I. The simulations assume 
spherical symmetry and use the cartoon code described 
in Appendix A. The unperturbed trumpet data is dis- 
cussed in Appendix B. The perturbations in the BSSN 
and gauge variables are defined by setting one of the in- 
coming modes to 



A. Incident mode Xe 



-0.0008 (r - ro)e-i2('--'-o)' 



(33) 



Here, 



X is one of xi, X2, Xs, X4 > 

X^, or Xe • Note that the vector and trace-free tensor 
modes 7 through 10 do not exist in spherical symmetry. 
Thus, the numerical tests performed here are restricted 
to modes 1 through 6. 

For the incoming modes X4 ; X5 Xe ' choose 
To = 10. For incoming modes Xij X2 and X3 we use ro = 
4.0. In each case we carefully examine the reflections in 
mode xt ■ This is the only mode in spherical symmetry 
that can propagate from the puncture boundary to the 
black hole exterior. 

Let me be specific about the definition of the charac- 
teristic fields in spherical symmetry. The definitions ( 24 ) 
include terms such as 9„5„„, which is shorthand notation 
for n'^daigbc'n^n'^)- The normal vector orthogonal to the 
r — const coordinate surfaces is n" = g""^ /^/(f^, which 
simplifies to S!^/ ^'g^r in spherical symmetry. Thus, we 
have 



Figure 
t = 3.4. ' 
ture r 



13 shows a perturbation in mode Xe at time 
^his pulse propagates inward toward the punc- 
^ .0, 



0. Figure 14 shows the incident pulse at t 
jus t before it hits the puncture boundary. Figures [T5| and 



161 show the reflected pulse Xe^ propagating outward at 
times 11.4 and 16.0, respectively. In each of these four 
figures, the dots (every fifth data point from a simulation 
with h — 1/50) lie on top of the solid curve (from a sim- 
ulation with h — 1/100). These results are convergent. 
In particular, the refiected wave pulse does not show any 
resolution-dependent time delay. 



T3 



T3 
O 



8e-05 
6e-05 
4e-05 
2e-05 h 


-2e-05 
-4e-05 
-6e-05 h 
-8e-05 
-0.0001 



time = 3.4 



1 



4 

r 








(34a) 


dn9AA 


= 2dr{gee / gee) / ^/g^ , 


(34b) 


A 




(34c) 




= {dr'f>)/^/g^ , 


(34d) 


K 


= K , 


(34e) 


A" 


= \/5^A'' , 


(34f) 


d„a 


= (9,.a)/v^ , 


(34g) 




= dr{y/g^P'^)/y/'g^ , 


(34h) 




= y/grrB'^ ■ 


(34i) 



where 9 is the usual polar angle in spherical coordinates. 
In the cartoon code described in Appendix A, the modes 
( 24 1 arc defined by using the frozen coefficients approx- 



imation to remove factors of the unperturbed fields grr 
and gee from the derivatives in Eqs. (34 1. For example, 
in Eqs. 

by {drgrrUig T 



and (25 1, the term 9„5„ 



is approximated 



The figures throughout this section use the following 
convention: The curves (solid or dashed, with various 
patterns) are obtained from simulations at resolution h = 
1/100. The data points (dots or crosses or other symbols) 
are obtained from simulations at resolution h — 1/50. 
Unless otherwise stated, only every fifth data point is 
displayed for the lower resolution case. 



FIG. 13: Incoming mode Xe, ^'t time 3.4. (The initial excita- 
tion is in Vfi.) 



0.0006 
0.0004 



(D 
T3 



to 

T3 
O 




FIG 
tion 



14: Incoming mode ^it time 8.0. (The initial excita- 
is in Xe •) 



It is informative to compare the amplitudes of the in- 
coming and reflected pulses. From Fig [13] we see that the 
incident mode Xe" has amplitude A ~ 8 x 10~^ at t = 3.4. 
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FIG. 15: Outgoing mode xt time 11.4. (The initial exci- 
tation is in Xe •) 
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FIG. 16: Outgoing mode xt time 16.0. (The initial exci- 
tation is in Xe ■) 



By t ~ 8.0 the amplitude has grown by almost an order 
of magnitude to A ^ 6 x 10~^. It continues to grow as 
the pulse approaches the puncture. The reflected mode 
^ initially has a very large amplitude. It drops rapid! 



to A 



4 X 10 at time t = 11.4, as shown in Fig. 



By t — 16.0 it has dropped an order of magnitude to 
A ^ 4 X 10~^. At equal distances from the puncture 
(say, r w 6.0), the reflected pulse is almost as large as 
the incident pulse. 

The small step in the data for xt seen in Fig 
r « 0.4 deserves some discussion. To begin, let us con 

xt 



at 



sider the modes Xa i xt xt before the pulse Xe ^i^s 



the puncture boundary. Figure [17] displays these modes 
at time t = 8.0. If our system of partial differential equa- 
tions were linear with constant wave speeds, these modes 
would not be excited at all. However, due to the nonlin- 
ear nature of Einstein's theory, these modes build in am- 



plitude as the initial pulse Xe propagates inward. Even 
the modes xt ' xt ^^'^ xt ' which are normally outgo- 
ing, develop nontrivial proflles that are carried inward 
along with Xe ■ The largest of these modes is xt ' ^'^^^ 
an amplitude A ~ 1 x 10~^. However, keep in mind 
that relative amplitudes among the different modes have 
little meaning since the normalization of the modes is 
arbitrary. 
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FIG. 17: Modes xji xt^ ^rid xt time 8.0, carried in- 
ward along with the initially excited mode Xe ■ Compare with 
Fig. [141 
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FIG. 18: Modes Xs > stud xt time 16.0. These waveforms 
are traveling inward. (The initial excitation is in Xe ■) 



The important feature of Fig. 17 is the small hump in 
X^ that peaks at about r = 3.8. This hump is initially 
part of the larger perturbation in X5 that is carried along 
with Xe ■ However, mode X5 naturally tends to move 
more slowly than Xe 1 ^ seen from the data in Table I. 
By t = 8.0 a small piece of X5 has broken free from the 
other perturbations and fallen behind. It continues to 
propagate inward, moving more and more slowly as it 



16 



approaches the puncture. Figure 18 shows the hump in 
at t = 16.0, just before it reaches the puncture. That 
figure includes a plot of xt ■ We see that the nonlinear 
couplings have created a new perturbation in Xe , that is 
carried inward along with Xs ■ This perturbation in xt 
is seen as a small step at r « 0.4 in both Figs. [T6| and [TS] 
Note that these features are all convergent. 

Finally, let us compare the incoming pulses Xe at < = 
3.4 and t — 8.0; these are pictured in Figs. 13 and 14 
Observe that the pulse e.t t — 8.0 is inverted relative to 
the pulse at t = 3.4. This is related to the breakdown 
in hyperbolicity that occurs when (1 — 2ae~'^'^) = 0. In 



particular, the factor (1 — 2ae ■*'^) appears in Eq. (24") 



for xt denominators of several terms in the 



inverse relations (25 1. This factor vanishes at r = 4.06 
for the trumpet geometry. More precisely, (1 — 2ae~*'^) 
is positive for r < 4.06 and negative for r > 4.06. Thus, 
as the pulse Xe passes through r — 4.06, the factors (1 — 
2ae-'^'^) in its definition (|24f) chang e sign. This change 
of sign is responsible for the inversion that is seen in the 
incoming pulse Xe ■ Figure 19 shows a time sequence of 



the inversion process. One can see from Figs. [T5|and|16| 
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FIG. 19: Time sequence showing mode Xe as it travels inward 
through the point r — 4.06 where hyperbohcity breaks down. 

that a similar inversion occurs in the reflected pulse xt 
as it travels outward and passes through r = 4.06. 



B. Incident mode X4 



Figures 20 



and 



21 



show the initial pulse X4 propa- 
gating inward at times t = 7.0 and t — 11.0, respec- 
tively. Figures 



22 



and 



23 



show a reflection in xt prop- 
agating outward~at times t — 11.0 and t ~ 15.0, respec- 
tively. All of these figures show nice convergence, with 
no resolution-dependent time delay. 

These results require further analysis an d explanation. 
Note that the pulse xt ' shown in Fig. 22 emerges from 
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FIG. 20: Incoming mode Xa a-t time 7.0. (The initial excita- 
tion is in xl ■) 
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FIG. 21: Incoming mode xl at time 11.0. The dots coincide 
with all of the data points for the low resolution simulation. 
(The initial excitation is in xl ■) 



the puncture boundary. Figure 21 shows that at time 
t = 11.0, the incident pulse is at r « 0.2. Evidently 
the perturbation in the outgoing mode xt is caused by 
something other than a reflection of Xi from the punc- 
ture boundary. 

we plot the modes xt^ xt ^ xt time 
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In Fig. 

t = 7.0. What is immediately clear from this flgure is that 
the nonlinear interactions have produced a large excita- 
tion in Xe ■ Observe that even with the "gauge system" 
described in Sec. III.D, the spatial variations in the char- 
acteristic speeds will introduce a source in the evolution 
equation for Xe that depends on Xa ■ This is a particu- 
larly strong coupling; note that the amplitude of Xe at 
t = 7.0 is several times larger than the amplitude of Xi- 
(Recall, however, that relative amplitudes depend on the 



the puncture region before the incident pulse X4 reaches normalization of the modes.) 
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FIG. 22: Outgoing mode xt time 11.0. (The initial exci- 
tation is in X4 •) 
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Compare with Fig. |20[ 
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FIG. 23: Outgoing mode xt at time 15.0. (The initial exci- 
tation is in X4 ■) 



Table I shows that near the puncture, mode Xe prop- 
agates more quickly than mode Xi- ^ consequence 
the excitation in Xq races ahead of X4 ■ One can see this 
already at time t = 7.0; Xe^ is centered at r w 1.9 while 
X4 is centered at r « 2.2. The large perturbation in Xe 
reaches the puncture boundary at t « 9.0 an d p rod uces 



22 



and 



23 



the large reflection in mode xt ' seen in Figs 

By time t — 11.0, the incident pulse X4 has reached 
r « 0.5. Observe that at t — 11.0, Fig. [22] shows a small 
hump in xt near r m 0.5. This hump is caused by the 
coupling between modes X4 ^rid Xe^- The tail of this 
hump is seen in Fig. [23] as an upturn in the data near the 
origin a,t t = 15.0. 
Figure 



25 



shows a close-up view of xt near r = at 
t — 15.0. This view reveals another feature, a bump in 
X^ that peaks around r w 0.6. Close inspection of the 
data shows that this bump has the following origin: The 



main reflection in xt creates an excitation in mode Xs^j 
which then propagates inward toward the origin. In turn, 
this incoming wave X5 creates a distortion in xt that is 
carried inward along with Xs^- It is this distortion that 
appears as the small bump in Fig. [25] 
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FIG. 25: Close-up vie w of outgoing mode xt at time 15.0. 
Compare with Fig. 23 (The initial excitation is in xl ■) 



The incident perturbation X4 reaches the puncture 
boundary at about t w 15.0. It does not produce a 
significant reflection in mode xt ■ Figure [26] shows the 
data for xt at t = 19.0, when such a reflection would 
be expected to appear in the range < r < 4.0. By 
this time the excitation in xt decayed to a "tail" 
near the boundary that oscillates as it flattens to zero. 



Note that the amplitude in Fig. 26 is quite small com- 
pared to the earlier reflections shown in Figs. ]22]andl23] 
More importantly, the results are convergent: there are 
no resolution-dependent features. 
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FIG. 26: Outgoing mode xt time 19.0. (The initial exci- 
tation is in X4 •) 



FIG. 28: Modes xf , X5 . and Xe" at time 11.0. The mode 
X^ is moving outward. The other modes are moving inward 
along with the initial perturbation in X5 ■ 



C. Incident mode Xs 



Figure 



27 



shows the mode Xs at time t = 11.0 as it 
propagates inward toward the puncture boundary. Non- 
hnear interactions cause the other modes to become ex- 
cited. By t = 7.0, mode Xe has been excited to an am- 
plitude of about A w 1 X 10^^. This mode moves more 
quickly than X5 and reaches the puncture at t ~ 9.0. 



about t « 20.0. It does not produce a distinct reflection 
in mode xt ■ Figure [29| shows mode xt at t = 22.0, when 
such a reflection would appear in the range < r < 4.0. 
The only excitation in xt at this time is a decaying tail, 
similar to the profile seen at late times when the initial 
perturbation is in X4 ; see Fig. 26 This tail is convergent. 



showing no resolution dependence. 
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FIG. 27: Incoming mode Xs at time 11.0. (The initial exci- 
tation is in X5 •) 



FIG. 29: Outgoing mode xt at time 22.0. (The initial exci- 
tation is in X5 ■) 
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at 



The modes Xi 7 xt 1 and xt are plotted in Fig 
time t = 11.0. The relatively large perturbation in xt ■ 
with amplitude A w 1 x 10~^, comes from the earlier re- 
flection of Xe ■ This mode is moving outward at t = 11.0, 
and is convergent. The other modes shown in Fig. [28] are 
moving inward, being carried by the initial perturbation 

The main pulse X5 reaches the puncture boundary at 



D. Incident modes xi) X2 and X3 

The results of these three cases are qualitatively very 
similar to one another. The incident pulse takes about 
t w 45 to propagate from its initial location at rg = 4.0 
to the puncture boundary. Long before that, the non- 
linear interactions create excitations in the other modes. 
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In particular, the excitation in mode Xe races ahead and 
reaches the puncture boundary at i « 4. As expected, 
this produces a reflection in mode xt that propagates 
outward to the black hole exterior. When the initial per- 
turbations in xii Xi or X3 tiit the puncture boundary, 
they do not appear to couple to any other modes. At 
late times we only see a residual tail in xX that decays 
in time. 



boundary, the mode develops a relatively large am- 
plitude profile. That profile decays in time, and does not 
appear to propagate outward. In Fig. 31 we see xX ^t 
time i ~ 55, after the initial pulse has hit the puncture 
boundary. This figure shows the residual profile and no 
resolution dependence. 



VI. INTERPRETATION AND DISCUSSION 
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FIG. 30: Mode xt ^t time 5.0. The three curves correspond 
to initial excitations in xi, Xa and X3- 
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Only the modes xt ^ind xt ahle to propagate from 
the puncture boundary to the black hole exterior. When 
starting the evolution from a wormhole configuration, 
modes xt ' xt' Xg ^^'^ xto outgoing at the punc- 
ture boundary. However, these modes quickly acquire 
negative speeds near the boundary and their character- 
istic curves never extend from the puncture boundary to 
the black hole exterior. 



The modes Xe f^nd Xs have constant speeds, it-^/g^ w 
±1 and ±y/3g^/2 « ±0.87 respectively, near the punc- 
ture boundary of a trumpet geometry. This implies that 
these modes propagate between the puncture boundary 
and the black hole exterior with no resolution depen- 
dence. These modes do not recognize that the puncture 
boundary moves when the resolution is changed. In ef- 
fect, modes xt '^^'^ xt propagate in the conformal geom- 



etry, not the physical geometry. 
The other modes, namely xi, X2, XSi 



xt 



xt 



XT, 



xt 



and Xu)' have speeds — Z?"" « — 0.53r near the puncture 
They propagate in the physical geometry. These modes 
are affected by the movement of the puncture boundary 
when the resolution is changed. Once the spatial slice 
settles to a trumpet geometry, all of these modes are 
incoming at the puncture boundary. 

The distinction between modes that propagate in the 
physical geometry and those that propagate in the con- 
formal geometry can be understood by examining the 
speeds from Tables I and II. The coordinate speeds for 
modes Xe^ and Xs^ built from the shift vector and 
the conformal metric component g^^ . The physical met- 
ric does not appear. On the other hand, the coordinate 

X^ and Xio are built from the 



± 



xt^ 



speeds for modes Xi 

shift vector and the physical metric component e^^'^g'"''. 

The time required for a perturbation in any mode to 
travel between the puncture boundary at r h/2 and a 
finite radius tq is 



FIG. 31: Mode xt time 55.0. The three curves correspond 
to initial excitations in xi, Xa and xs- 



Figure 
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I shows Xe for the three cases with initial exci- 
tations in xi, X2 and X3 at time t = 5.0. The large pulses 
near r — 1.2 are reflections from x^ ■ These wave forms 
are traveling outward with positive speed. The smaller 
pulses near r = 3.6 are excitations that are carried in- 
ward along with the initial perturbations in xi, X2 or xs- 
None of these curves show any resolution dependence. 
As the initial perturbation passes through the puncture 



T: 



dr 

h/2 \v\ 



(35) 



where v is the coordinate speed. For the modes that 
propagate in the conformal geometry, v is constant near 
the puncture boundary and T depends linearly on the 
grid resolution h. In the limit /i — > 0, this dependence 
drops out. In this sense, the propagation time does not 
depend on resolution. For the modes that propagate in 
the physical geometry, v behaves like |t;| ^ r near the 
puncture boundary of a trumpet. For these modes T 
diverges like |ln(/i)|. As the resolution is increased h is 
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decreased and the puncture boundary is pushed farther 
down the trumpet throat. The propagation time T goes 
to infinity. 

Once the wormhole has evolved into a trumpet, all of 
the modes that propagate in the physical geometry have 
negative speeds at the puncture boundary. The concern 
is that a disturbance in one of these modes could prop- 
agate inward to the puncture boundary where it might 
couple to one of the outgoing superluminal modes xt "-"^ 
■ Since the time required for the initial disturbance 
to reach the puncture boundary is resolution dependent, 
the reflection would be delayed in time as the resolution 
is increased. 

The numerical experiments in Sec. V suggest that there 
is no coupling at the puncture boundary between the 
modes that propagate in the physical geometry and the 
modes that propagate in the conformal geometry. These 
experiments were limited to spherical symmetry, so not 
all possible couplings could be tested. 

The lack of coupling between the two types of modes 
can be understood by considering the finite differencing 
stencil as it appears from the perspectives of the con- 
formal and physical geometries. From the point of view 
of the physical geometry, the finite differencing stencil is 
one sided at the puncture boundary. We can see this in 
Fig. [2j For any point on the puncture boundary, the legs 
of a stencil will always extend outward (into the com- 
putational domain). For the modes that travel in the 
physical geometry, the puncture method correctly avoids 
placing any boundary conditions at the puncture bound- 
ary by using one-sided stencils. These modes are sim- 
ply advected through the puncture boundary and off the 
computational domain. 

For the modes that travel in the conformal geometry, 
the finite differencing stencils appear as standard stencils 
that surround the origin of a computational grid in R^. 
For these modes, there is no boundary and the puncture 
method correctly treats the puncture like any other point 
in the computational domain. 

Because the modes 'X^ propagate in the confor- 

mal geometry, which is smooth at the origin, we expect 
the puncture method to impose some type of smoothness 
conditions at r = 0. It is not entirely clear what those 
conditions might be since the unperturbed fields are not 
all smooth at that point. The following observations are 
suggestive. If we assume that the perturbations in the 
vector fields i?" and Z?** are smooth at the orig in, then 
we must have = and /3" = 0. From Eq. (24 ") we see 

that the modes xt obey Xe^fi) = TdnP^i^) at r = for 
a trumpet geometry. Thus, we expect that the relation 



X^(0)+X6"(0)-0 



(36) 



should hold for all of the simulations in Sec. V. A close 
inspection of the data shows that this is indeed the case. 

For BSSN with standard gauge conditions, some modes 
effectively propagate in the physical geometry and some 
effectively propagate in the conformal geometry. In light 



of this understanding, let us examine the scalar field ex- 
ample of Sec. II. In that case the characteristic speeds 
are v — near the puncture boundary. 

Both the incoming and outgoing modes propagate in the 
physical geometry. The time T for these modes to travel 
between the puncture boundary and any finite location 
diverges like | ln(/i)|. Moreover, the finite difference sten- 
cil for both of these modes is one sided at the puncture 
boundary. This is bad. For a properly constructed nu- 
merical code, the outgoing mode (the mode that is prop- 
agating into the computational domain) should be fixed 
by boundary conditions. This does not happen when we 
naively apply the puncture method as in Sec. II. Corre- 
spondingly, the code used in Sec. II does not work — the 
simulations do not converge to a continuum solution as 
the resolution is increased. 



On the other hand, a properly constructed numerical 
code for BSSN with standard gauge should have precisely 
the properties of the puncture method. That is, the code 
should use something like one-sided finite difference sten- 
cils for the modes that propagate in the physical geome- 
try, since all of these modes are incoming at the puncture 
boundary. For the modes that propagate in the confor- 
mal geometry there is no boundary, so the code should 
use the same stencils near the origin as elsewhere. 

The analysis in this paper assumes that the gamma- 
driver shift condition includes all advection terms. There 
are three such terms, one in Eq. ( 10 1) and two in 



Eqs. ( 10 i). Let me comment on the cases in which one or 



more advection terms are omitted. What we look for is a 
case with an outgoing mode whose coordinate speed goes 
to zero at the puncture. To determine the speeds near 
the puncture boundary, I will assume that the geometry 
evolves to a stationary trumpet with a ~ ar, (5^ « hr, 
g^"^ ~ g and ip — In(pr) /2 near the puncture boundary, 
where a, b, g and p are all constants. This assumption 
holds if we include all advection terms, and it holds if we 
include no advection terms. 



With the assumption above, there is only one scalar 
mode and one vector mode whose speeds can approach 
zero through positive values near r — 0. This happens 
when the second advection term in Eq. i\10\ i) is dropped 
but the other advection terms are kept. [The second 
advection term is contained in (9tA°)i.hs-] In that case 
the scalar mode has speed b^r^ /g and the vector mode 
has speed Ab'^r^ /{'3g) near the puncture boundary. If ei- 
ther of these modes is excited at the puncture bound- 
ary, it will propagate to the black hole exterior with 
a resolution-dependent time delay. More precisely, the 
propagation time will diverge like 1 //i^ as the resolution 
is increased. The conclusion is that we should not drop 
the second advection term in the gamma-driver shift 
equation Eq. (lOi) unless one or both of the other ad- 



vection terms are dropped as well. 



21 



APPENDIX A: CARTOON BSSN CODE 

The cartoon code evolves the BSSN plus standard 
gauge system in spherical symmetry using Cartesian fi- 
nite difference stencils. The computational grid is shown 
in Fig. |32] The puncture is at the origin. The phys- 
ical grid consists of the single line of grid points at 
x/h = y/h = 1/2 and z/y = 1/2,3/2,..., where h is 
the grid spacing. These points are shown as solid dots 
in the figure. The open circles represent buffer points, 



o : o : o o : o 



o o o o o 



FIG. 32: Grid for the cartoon code. The physical grid points 
are filled dots, and buffer points are open circles. 

which are needed for computing finite difference deriva- 
tives. The figure shows one layer of buffer points; the 
actual code uses two layers of buffer points. 

The buffer points are filled by spherical symmetry us- 
ing known values of the fields at physical grid points. 
Consider a scalar field such as the lapse function a. 
The value of a at a buffer point {Xg,yg,Zg) is found 
by considering the coordinate distance from the origin: 

rg = ^Xg + yg + z'^. We compute the value of z along 
the physical grid line for the point at radius Vg. That is, 
we set {h/2Y + {.h/2Y + z"^ — Vg and solve for z\ this 
gives 



xl + yl + zl 



h/2 



(Al) 



I use an eighth-order interpolation along the z direction 
to find a at radius rg. 

As an example, consider the buffer point 
(5/i/2,3ft./2,9/i/2). According to Eq. (|AT) we have 



z = 5.315/i. Thus, the point {h/2, h/2, 5.315/i) lies along 
the physical grid line at the same radius as the buffer 
point. Note that z — 5.315ft. lies between 9/i/2 and 
llft/2. The value of a at {h/2,h/2,b.?,lbh) is found 
by interpolation over the eight physical points with 
z/h = 3/2, . . . 17/2. This is the value assigned to a at 
the buffer point. 

The conformal connection A°, the shift P"" , and the 
auxiliary field B"" are all spatial vectors. Let V"" denote 
one of these vectors. In spherical symmetry, the Carte- 
sian components of any vector field can be written as 



y° = s{ry- 



(A2) 



where S{r) is a scalar field. Here, x"" are C arte sian co- 
ordinates and r — x"^ + y"^ + z"^ . Relation ( A2 ) can be 
inverted. 



S=V'' 



(A3) 



where indices on x"" are lowered with the flat Cartesian 
metric 5ab- In filling buffer points for vectors, we first 
apply Eq. (A3) to compute the scalar S. Next, we apply 



the scheme outlined above to fill the buffer points for 
S. Finally, the buffer points for V"" are computed from 



Eq. (A2l 



Let Tab denote one of the symmetric tensor fields gab 
or Aab- In spherical symmetry we have 



Tab^P{rf-^ + Q{r)5ab 



(A4) 



where P and Q are scalars. This relation can also be 
expressed as 

Tabdx^dx^ = {P + Q)dr^ + r^QdQ^ , (A5) 

where dfi^ is the metric on the unit sphere. Buffer points 
for Tab are filled by first computing the scalars from the 
inverse relations 



P + 3Q = TabS"'' 



P+Q = Ta, 



x'^x^ 



(A6a) 
(A6b) 



We then fill the buffer points for the scalars, as outlined 
above. Finally we apply Eq. (A4| to obtain buffer point 
values for Tab- 

Note that the calculations discussed here rely on the 
fact that all of the fields transform as simple tensors with 
no density weights. This is the case for the covariant 
formulation of BSSN and the standard gauge P2]. 

The cartoon BSSN code uses fourth-order Runge- 
Kutta for time integration. It uses standard fourth-order 
centered finite difference stencils for spatial derivatives, 
like those displayed in Eqs. (|6|. The only exceptions to 
this rule are the advection terms, which have the form 
/3°'daF for a tensor or tensor component F. These terms 
are upwinded. The code as it is now written includes two 
layers of buffer points. This limits the size of the finite 
difference stencils such that the advection terms are only 
third-order accurate. For example, the finite difference 
approximation to j3^dzF with > is 



(/3^a,F)fc = ^/3,^(-2Ffe_i 



^3Ffc + 6Ffc+i-Ffe+2) (A7) 



where k labels the grid points and h is the grid spacing. 
The overall accuracy of the cartoon code is limited to 
third order by the presence of these advection terms. 



Figures 33 and 34 show the results of convergence tests 
for the cartoon code. The data for both of these fig- 
ures come from a simulation in which a wormhole slice 
of Schwarzschild evolves into a trumpet geometry. The 
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data are taken at time t = 10. At this time the various 
fields, including the lapse function and shift vector, have 
settled fairly close to their trumpet values in the imme- 
diate vicinity of the puncture (say, r<0.1) but are still 
changing rapidly away from the puncture (r>0.1). 

Simulations at five different grid resolutions were used 
for these tests. Let us label these resolutions by integer 
subscripts 1, 2, 4, 8, and 16. The grid spacings are hi ^ 
1/12.5, h2 = 1/25, /i4 = 1/50, hs = 1/100, and hie = 
1/200. 

Figure [33] shows the absolute value of the Hamiltonian 
constraint. The curves in this figure are scaled by factors 
of 8 at successive resolutions, appropriate for a third- 
order numerical scheme. As one can see, the curves over- 
lap in the limit of increasing resolution. This shows that 
the code is indeed third-order accurate. 
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FIG. 33: Absolute value of the Hamiltonian constraint at five 
different resolutions. The curves are scaled by powers of 8 to 
show third-order convergence. 



FIG. 34: Absolute value of the difference between the confor- 
rnal factor at successive resolutions. The curves are scaled by 
powers of 8 to show third-order convergence. 



three higher resolution curves overlap nicely. Thus, the 
code appears to be third-order convergent at r = 0.3. 

Similarly, let us consider the point r = 0.1. None of the 
curves overlap at this point. This is due to insufficient 
resolution. It is fairly clear from examining the sequence 
of curves that if we were to add a higher resolution run, 
with = 1/400, that the curves 841^161 and 8^\n32\ 
would closely agree. We then expect that at r = 0.1 the 
code is third-order convergent. In fact, the code appears 
to be third-order convergent at all points r > 0. This 
is as it should be since the data are smooth everywhere 
except at r = 0. A properly constructed code should 
be convergent at all points in the computational domain, 
excluding only the puncture points. 



Figure[34]shows the result of a three-point convergence 
test for the conformal factor (p. Each curve is obtained 
from the difference between the conformal factor at suc- 
cessive resolutions, scaled by an appropriate power of 8. 
The curves in this plot also overlap in the limit of increas- 
ing resolution. This confirms that the code is third-order 
convergent. 

These results require some elaboration. It has been 
stated in a number of publications that puncture evolu- 
tion codes do not converge near the puncture. What we 
see from these figures is that the errors are large near 
the puncture. This does not imply a lack of convergence. 
A code is convergent if, at any given point in the com- 
putational domain, the scaled errors coincide with one 
another in the limit of high resolution. For example, con- 
sider the point r = 0.3 on the graph of Fig.[33j This point 
is about half way between 0.1 and 1 on the logarithmic 
scale. The two lowest resolution curves do not coincide 
with the higher resolution curves at r = 0.3, due to the 
presence of large finite differencing errors. However, the 



APPENDIX B: TRUMPET GEOMETRY IN 
GAMMA-DRIVER COORDINATES 



In Sec. V we considered the evolution of perturbations 
on a trumpet slice of a single Schwarzschild black hole. 
The perturbations were defined in terms of characteristic 
fields through Eqs. (25). This appendix is devoted to a 



description of the unperturbed trumpet geometry. 

In this paper, the term "trumpet geometry" refers to 
a stationary 1+log slice of the Schwarzschild geometry. 
Such a slice can be expressed in various spatial coordi- 
nate systems. One technique for finding a trumpet slice 
in isotropic coordinates has been described elsewhere 
[lOl 124] , What is needed for the analysis in Sec. V is 
a trumpet slice in stationary gamma-driver coordinates. 
That is, the data should satisfy the gamma-driver shift 
equations ( 10 i,h) including advection terms and with 
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time derivatives set to zero: 



= P^DcP"" + -B" , 







(Bla) 
(Bib) 



The damping parameter rj was set to zero in all simu- 
lations presented in this paper; it is set to zero here as 
well. 

Let me begin by reviewing the construction of trumpet 
data in isotropic coordinates, then show how the results 
can be extended to stationary gamma-driver coordinates. 
A stationary, spherically symmetric geometry can be de- 
scribed by the metric 

ds^ = ~{a^ + (3^)dt^ + 2/3 d£ dt + df + R'^dn'^ (B2) 

where the lapse a, shift /3, and areal radius R are func- 
tions of the proper distance coordinate £. Recall that 
dfl^ = dO^ + sin^ 6 d<p is the line element for the unit 
sphere. The Einstein equations imply 



(3^ = \- 2M/R , 



(B3) 



where M is an integration constant and the prime de- 
notes d/d£. The 1-Klog slicing condition ( 10 '), along with 
stationarity, gives 



2/r 



AR' 



(B4) 



Equations ( B3 1 can be used to eliminate a and /3 from 



Eq. (B4l; this yields 



R" 



2R'[3M -2R + 2RR'^] 
R[2M - R- 2RR' + RR'^] 



(B5) 



In Ref. [9], Hannam et al. argue that the numerator and 
denominator on the right-hand side of Eq. ( B5 ) must 



vanish separately at some particular value of I. Solving 
these conditions simultaneously shows that the equations 
R' = -2, + yiO and M/R = 4(-3 -I- \/T0) must hold for 
some value of t. 



Equation (B4) can be integrated to /3i?^ — ■\/C'e"/^, 



where C is an integration constant. By using Eqs. (B3) 
we can rewrite this expression in terms of R and R' : 



Ce^' ^ R^iR'" -1 + 2M/R] 



54rD'2 



(B6) 



The values previously obtained for R and R' at one value 
of £ can be used to solve for the constant, with the result 



C 



4„3- 



/lO 



128(VT0- 3)s 



1.55 



(B7) 



Note that in the limit as £ —oo, the trumpet geometry 
satisfies R' 0. Equation (B6) has two real solutions 



for R in this limit, but only one has R' approaching 
through positive values. For that solution we find that 
the trumpet throat has radius R w 1.312 Af in the limit 
f ^ -oo. 



The proper distance can be written as 

d^ 
' R' ' 



(B8) 



where R' is considered to be a function of R. This de- 
pendence, i?' as a function of R, is found by applying 
Newton's method to Eq. (|B6|). Then the integral (|B8| 



can be evaluated numerically, starting from large nega- 
tive £ where R « 1.312 Af. This gives £ as a function of 
R. 

Equivalently, the analysis above defines i? as a function 
of £. The data a and as functions of £, are found from 
Eqs. (B3). The components of the extrinsic curvature 



are computed from the spacetime metric ( B2 1 as 



Km 



Rf^_2^ 
2R'~R 

I3R , K^^ = l3Rsin^0 



(B9a) 
(B9b) 



These are determined as functions of proper distance £ 
from the numerical solution R{£). 

We now switch to isotropic coordinates, with spatial 
metric 



ds^ = ^"^{dp^ + p^dn^) 



(BIO) 



Comparing with the spatial part of the metric ( B2 1 , we 
have ^ = and dp/di = p/R. Since R is known 

as a function of £, this second relation can be integrated 
numerically to give p as a function of t. Turning this 
around, we can consider £ ss a, function of p. Then 
the conformal factor is determined as a function of 
the isotropic coordinate p: expHcitly, ^' — yj R{£{p))/ p. 
Likewise, the lapse function becomes a known function 
of the isotropic radius: a — a{£{py). The shift vec- 
tor in isotropic coordinates is 13^ = {p/R) (3, and the 
radial-radial component of the extrinsic curvature is 
Kpp = {R/ pYKu. These are now determined as func- 
tions of p since the dependence £{p) is known. 

The stationary gamma-driver shift equations with rj = 
are written in Eqs. (Bll above. Assuming 0, the 



second of these equations can be integrated to obtain 
j^a = _|_ const. We can choose initial data such that 
= A"". Then the first of Eqs. (|B1]) and the definition 
^ for A° yield 







r.bC 



be 



(Bll) 



Our goal is to find a change of spatial coordinates, from 
isotropic coordinates p, 9, (f) to "gamma-driver coordi- 
nates" r, 9, (p, such that Eq. (Bll I is satisfied in the new 
coordinate system. 



We begin by rewriting the spatial metric (BIO) as 



ds' = 



2/3 



4/3 

^ dr' 



2 \ 2/3 



dn' 



(B12) 
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Here, p' is the derivative of p with respect to the new 
radial coordinate r. Observe that the conformal part 
of this metric (the factor in square brackets) has deter- 
minant r^s\T?Q. With this choice the conformal metric 
will have determinant 1 in Cartesian coordinates. Now 
use the conformal metric from Eq. (B12l to write the 



gamma-driver shift equation (Bill explicitly 



1 

3^ 



rp 
P 



2/3 



4p 
rp' 



Pn" 



2p' ip'^p" 



(B13) 



Here, the shift has functional dependence = p''{p{r)). 



Equation (B13l is a second order differential equation 



for /9 as a function of r. We can turn this around and view 
it as an equation for r as a function of p. In doing so, 
we must make the changes p' = 1/r' and p" = —r"/p'^, 
where primes on r denote derivatives with respect to p. 
With this view of Eq. (B13l, the shift has functional de- 



pendence P'' — P''{p)- Recall that P^ is known as a func- 
tion of p from the numerical trumpet solution in isotropic 
coordinates. 

The analysis above yields the following second order 
elliptic differential equation for r as a function of p: 



ipr'y^' + 2{pr'fl\ - 'i{pr'fl\'' - 2ppPpP'r'\"^ 

I 



2ppP\'r"'' -2{p\'f'\ 



(B14) 



It is helpful to make the change of variables a = r/ p. The 
resulting equation can be solved for cr as a function of p 
with boundary conditions cr = 1 as p — > and p — > cx). 
Any number of techniques can be used; I use a simple 
multigrid relaxation scheme. Once we have found r as a 
function of p, the usual tensor transformation equations 
can be used to determine the lapse function, shift vector, 
and extrinsic curvature in the new coordinates r, 0, and 

I have written a stand-alone numerical code to carry 
out the above calculations. Figure |35]is a plot of a versus 
p. It is immediately clear that a remains close to unity. 
Consequently the stationary gamma-driver coordinate r 
differs by a relatively small amount from the isotropic 
radius p. 
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FIG. 35: The ratio a = r/p, where r is the radial coordinate 
for the gamma-driver shift coordinates and p is the isotropic 
radius. 



nents of the conformal metric in the stationary gamma- 
driver coordinates. It is difficult to tell from this figure, 
but both components have a vanishing slope at r = 0. 
Figure [37| is a plot of the lapse function and Fig.[38|shows 
the radial components of the shift vector and the auxil- 
iary field. 



c 

0) 

c 
o 

Q. 

E 
o 
o 
o 

■•— » 
03 




Figure 36 is a plot of the radial and angular compo- 



FIG. 36: Conformal metric components grr and gee/r"^ for 
the trumpet geometry in gamma-driver shift coordinates. 

The coordinate location of the black hole horizon in 
stationary gamma-driver shift coordinates is r = 0.842. 
(Again, as in the body of this paper, I have set M — 1.) 
In isotropic coordinates, the horizon location is p = 
0.839. The location where e^''' — 2a and hyperbolic- 
ity breaks down is r = 4.06 in stationary gamma-driver 
coordinates. To three significant digits, this location has 
the same value in isotropic coordinates: p = 4.06. 

My original strategy for this research project was to use 
the stand-alone code to generate a very high resolution 
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FIG. 37: Lapse function a for the trumpet geometry in 
gamma-driver shift coordinates. 
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FIG. 38: Shift vector and auxiliary field components (3^ and 
for the trumpet geometry in gamma-driver shift coordi- 
nates. 



data set for the unperturbed trumpet. That data would 
then be read into the cartoon BSSN code and interpo- 
lated onto the grid used for evolution. This strategy does 
not work very well. One of the difficulties is that there 
are many places in which numerical errors can accumu- 
late. In particular, we have the integration in Eq. ( B8 1 , 



the integra tion of dp/di = p/R, the solution of the el- 
liptic Eq. (B14l, the transformations from isotropic to 



stationary gamma-driver coordinates, and finally the in- 
terpolations onto the evolution grid. With my particular 
implementation, the errors in the unperturbed trumpet 
data were large enough to keep the evolution code from 
showing clean third-order convergence. 

Another way that one can generate a trumpet geome- 
try in gamma-driver coordinates is to begin with worm- 
hole data and evolve that data until it settles to a sta- 
tionary state. When a wormhole evolves to a trumpet. 



there is an adjustment pulse that propagates out from the 
puncture boundary to infinity. It takes a time of about 
t = 100 for the trailing edge of this pulse to propagate 
beyond r = 50. To be sure that the data are not contam- 
inated with errors from the outer boundary, I evolve the 
initial wormhole for a time oi t = 110 on a grid with an 
outer boundary at r^^ax = 170. The data beyond r = 50 
is then discarded. 

The unperturbed data produced by this second 
method, evolving a wormhole for t — 110, works fairly 
well for studies with initial perturbations in modes X4 ; 
X^, and ■ In these cases, with tq = 10, the pertur- 
bations reach the puncture boundary within a time of 
i « 20 or less. A total simulation time of around t — 25 
is sufficient to reveal any reflections that might arise in 
mode xt ■ With the outer boundary at r^ax = 50, we 
can be confident that the results are not contaminated 
with errors from the outer boundary. 

The studies with initial perturbations in xii X2i Sind 
X3 require much longer run times, up to i = 120. This is 
because these modes move very slowly toward the punc- 
ture. For these cases, the second method of generating 
unperturbed trumpet data would require us to evolve 
the wormhole for t « 180 on a grid that extends to 
'''max ~ 300. This is not feasible with my current code. 
Thus, in Sec. V.D I used the first method, the method of 



Eqs. (Bl) through (B14), to generate unperturbed trum- 
pet data. 

The two methods of calculating unperturbed trumpet 
data yield similar results. In fact, one cannot tell the dif- 
ference between the two methods by examining graphs 
such as those in Figs. 36 — 38 The differences only ap- 



pear upon more careful examination, such as convergence 
testing. 

Another problem with the first method is that, near 
the puncture, the data is in some sense "too good". For 
much of my work I used the first method to generate 
trumpet data with a resolution of Ar w 0.00069. That 
data was interpolated onto an evolution grid with a typi- 
cal resolution of around h = 0.01. With the cartoon evo- 
lution code, as with any puncture method code, the res- 
olution is very poor near the puncture boundary and the 
truncation errors are high. On the other hand, the trun- 
cation errors in the data produced by the first method 
are relatively low near the puncture boundary. Because 
of this mismatch, the unperturbed trumpet data is not 
quite stationary when it is evolved with the evolution 
code. In particular, the grid points near the puncture 
boundary evolve as they adjust to the larger truncation 
errors of the evolution code. This adjustment creates a 
perturbation that propagates outward through the com- 
putational domain. The adjustment is difficult to notice 
on a plot of the BSSN or gauge variables. However, the 
adjustment readily appears as an excitation in the char- 
acteristic fields near the puncture. In some cases the 
excitation is much larger than the perturbations we wish 
to study. 

I can correct for this problem in the following way. For 
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each of the simulations with initial perturbations in xi, 
X2 and X3, I "normalize" the data by subtracting the 
results obtained from a simulation with no perturbation. 
This removes the adjustment pulse from the data. The 
normalized data is easier to interpret than the raw data. 

The same normalization scheme is applied to the simu- 
lations with initial perturbations in , X5 and Xe > but 
for a somewhat different reason. In these cases the unper- 
turbed trumpet data is derived using the second method, 
the method of evolving a wormhole for t = 110. However, 
the mode analysis is very sensitive to any change in the 
stationary trumpet data. What can be seen in the mode 
plots is that the unperturbed trumpet data is not en- 
tirely stationary. Even after a run time of t = 110 there 
is still a small amount of evolution taking place near the 
puncture boundary. For example, when the unperturbed 
trumpet data (prepared by the second method) is evolved 
for another t = 20, the value of mode Xe drifts from 
to —2.5 X 10~^ near the puncture boundary. These ef- 
fects are removed from the data for perturbed trumpet 



simulations by subtracting the data for an unperturbed 
trumpet simulation. 

To summarize, the simulations from Sec. V that have 

initial perturbations in modes X4 i X5 a-^d Xe ^^'^ ^^"^ 
second method for generating unperturbed trumpet data. 
The simulations that have initial perturbations in modes 
Xi, X2 and X3 use the first method for generating un- 
perturbed trumpet data. In all cases, the data shown 
in Sec. V have been normalized by subtracting the data 
obtained from a simulation of the unperturbed trumpet. 
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